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Quantum-ionic features in the absorption spectra of homonuclear diatomic molecules 


A. Crawford-UrangaQ D. J. MowbrayQ and D. M. CartlamoncH 

Nano-Bio Spectroscopy group and ETSF Scientific Development Center, 

Departamento de Fisica de Materiales, Centro de Fisica de Materiales CSIC-MPC and DIPC, 
Universidad del Pais Vasco UPV/EHU, Avenida de Tolosa 72, E-20018, San Sebastian, Spain 

(Dated: February 3, 2015) 

We show that additional features can emerge in the linear absorption spectra of homonuclear diatomic 
molecules when the ions are described quantum mechanically. In particular, the widths and energies of the 
peaks in the optical spectra change with the initial configuration, mass, and charge of the molecule. We intro¬ 
duce a model that can describe these features and we provide a quantitative analysis of the resulting peak energy 
shifts and width broadenings as a function of the mass. 
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I. INTRODUCTION 

Molecular spectroscopy deals with the response of a 
molecule interacting with an external electromagnetic field. 
The development of attosecond sources M allows one to 
probe in real time coupled electron-ion dynamics after pho¬ 
toionization processes. Two types of processes are seen in 
these experiments. The motion of the ions is associated with 
chemical transformations such as dissociation 0] in the fem¬ 
tosecond domain. The motion of the electrons is associated 
with electronic rearrangement processes such as charge redis¬ 
tribution Si, localization 13, i as well as ionization pro¬ 
cesses such as tunneling Si in the attosecond domain. 

Modeling coupled electronic-ionic dynamics in photoion¬ 
ization processes is a formidable challenge for most systems. 
For this reason, previous studies have been limited to one (Ht) 
and two (FF) electron benchmark systems SdMi. a fun 
coupled electronic-ionic 3D treatment has only been achieved 
for the one electron system FI), where the ionic motion is con¬ 
fined to the direction of the laser’s polarization GUB- For 
a full quantum mechanical treatment of two electron two ion 
systems (FF), it is necessary to confine both the electronic 
and ionic motion to the laser’s polarization direction. This 
is a reasonable semiclassical approximation, as the electronic 
and ionic motion should be predominantly along this direction 
[0. Therefore, for most molecules, any quantum-ionic fea¬ 
tures are typically neglected by instead using classical approx¬ 
imations, e.g., the Born Oppenheimer approximation (BOA) 
and Ehrenfest dynamics (ED). These approaches rely on a 
weak coupling between the electronic and ionic wave func¬ 
tions. However, the validity of such approximations breaks 
down for light atoms, when hybridization between the elec¬ 
tronic and ionic wave functions must be included. A quan¬ 
tum versus classical treatment of the ions has been previously 
used to investigate the localization JUl, nonsequential double 
ionization in and harmonic generation |Qj]] of FF, as well as 
the dissociation [J2] and proton kinetic energies |13] for H(. 
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The aim of this paper is a comparison between a quantum 
mechanical (QMI) and classical (BOA/ED) treatment of the 
ionic motion to describe coupled electronic and ionic pro¬ 
cesses JH. In particular, we consider three and four body 
systems of electrons and ions for which a fully quantum me¬ 
chanical treatment of the coupled electron-ion system is fea¬ 
sible. This comparison with respect to the QMI solution is 
performed both for the static spectra and for the time depen¬ 
dent linear response spectra. In fact, we find significant differ¬ 
ences between the QMI and BOA/ED spectra. These features 
can be quantitatively analyzed using a simple two-level two- 
parameter model based on the BOA electronic energy levels 
and the electron-ion mass ratio. The results of our work will 
help us to determine the domain of applicability of the simpli¬ 
fied BOA and ED approaches to interpret coupled electron-ion 
experiments for more complicated systems. 

The paper is organized as follows: in Sec. [II] we introduce 
the theoretical methods and models employed to simulate the 
coupled electronic and ionic processes; in Sec. m we explain 
the methodology to obtain both the ground state and time de¬ 
pendent linear response spectra, as well as the computational 
details of our calculations; in Sec.[IV] we show our results for 
both, the Hj and H 2 molecules, which we then analyze ac¬ 
cording to the model we provide; and finally, in Sec. [V] we 
summarize the main conclusions and relevant results of our 
work. Atomic units a.u. (h = m e = e = ao = 1) are used 
throughout, unless stated otherwise. 


II. THEORETICAL BACKGROUND 

A. Quantum electron-ion approach 

A many-body system composed of N ions and n electrons, 
where both the electrons and ions are treated quantum me¬ 
chanically (QMI), is described by the total electron-ion time- 
dependent Hamiltonian 

H(t) = fj + f e + V„ + V Ie + V ee + V ext (f), (1) 

where 7/ and T e are the ionic and electronic kinetic energy 
operators, respectively, and V u , V te , V ee , and V ext (t) are the 
ion-ion, ion-electron, electron-electron, and external potential 
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energy operators, respectively. The kinetic energy operators 
take the form 


N , 


-V 2 

v or’ 


( 2 ) 


where A7„ is the mass of ion a, and 

f-t-H 


i= 1 


The interaction between the ions is given by 


v n = - V 


a±p 


(3) 


(4) 


where (2/3: R<> and R /; are the corresponding charges and 
positions of ion a and [i. Similarly, the electron-electron re¬ 
pulsion is 


V,, 



ij= 1 
i*j 



(5) 


where r, and r ; are the positions of electrons i and j, while the 
interaction between electrons and ions is 


N n 


a— 1 i— 1 


Qa 

|r,- - Rffl' 


( 6 ) 


Finally, V ext (t) describes the interaction of the system of 
electrons and ions with an external electromagnetic time- 
dependent field, defined explicitly in Sec. lIIIBl 

The time-dependent Schrodinger equation takes the form 

dtf/ d 

i— = i—it/(RiSi,R 2 S 2 , ...R N SN',fisi,r 2 s 2 ...r n s n ,t) 
at at (7) 

= H(t)i//(RxS i, RiS 2 , -R nS jv; r i H , r 2 s 2 ...r n s n , t). 


where <// is the time-dependent electron-ion wavefunction. 
This depends on the positions R, and r, and on the spin coor¬ 
dinates S a and Si of ion a and electron i, respectively. 

For time-independent problems ( V ext (t ) = 0), the general 
solution of the time-dependent Schrodinger equation can be 
written as 

1/3 = ^ Cke~ ,ett i//k(RiS i, R2S 2, ...RjvS n\ ^isi,r 2 s 2 ...r n s n ) 

k 

( 8 ) 

where Sk and il/k are the k' h eigenvalue and eigenstate of the 
electron-ion stationary Schrodinger equation 

Hi//k = Hi// k (R\S i,R 2 5 2 , ...R N S N -,risur 2 s 2 ...r„s n ) 

= Sklf/k(RlS 1 .R 2 S 2 , ...Rjv5/v;riSl,r2S2—rnSn), 


H = f, + f e + Vn + V Ie + V ee . (10) 


We will next focus on the time-independent solution until in¬ 
troducing an external field in Sec. lIIIBl 

Solving the QMI problem is very demanding computation¬ 
ally for many-body systems. In fact, it quickly becomes un¬ 
feasible for systems with more than three independent vari¬ 
ables. For this reason, we restrict consideration herein to one 
or two-electron diatomic molecules whose motion is confined 
to one direction (see Sec. Ill Dk By applying an appropriate co¬ 
ordinate transformation, such systems may be modeled with 
only two or three independent variables (see Appendix^. In 
Secs. lllBl and lHCl we introduce two of the most widely used 
approximations to simplify the general many-body electron- 
ion problem. 


B. Born-Oppenheimer approximation 


Within the Born-Oppenheimer approximation (BOA) [16|, 
the total electronic-ionic wavefunction i/r is assumed to be sep¬ 
arable into an ionic x an d electronic ip part. As the electrons 
move much faster than the ions, we assume that the kinetic en¬ 
ergy of the ions does not cause the excitation of the electrons 
to another electronic state, i.e., an adiabatic approximation. 
Such an approximation is valid so long as the ratio of vibra¬ 
tional to electronic energies, E v n, to E e i ec , which goes as the 
root of the electron-ion mass ratio, i.e., E v u,/E e \ ec « yfmJM, 
is small |16|(see Appendix iBl for details). Since for a proton 
M p ~ 1836m f and E vih /E e j ec ~ 0.02, the BOA is expected 
to work quite well for our molecules. We thus may neglect 
77 = 0 from Eq. (I Kil l, although the electrons still feel the static 
field of the ions ( V e i , Vn + 0). 

The separable BOA solution 1 ft of the electron-ion stationary 
Schrodinger equation <[8]) is given by |12] 

=T( R i‘Si: R 252,-RA35/vV Rl ’ R2 ’' R ' v) (riSi,r 2 S2-r„s„), 

( 11 ) 

where x depends on the ionic coordinates only and p depends 
on both the electronic coordinates and on the ionic coordinates 
which, however, only enter into the electronic wavefunctions 
as parameters. As shown in Ref. El, this may be done for 
the Hamiltonian of Eq.[[| without loss of generality. 

If we insert Eq. (fill directly into Eq. <[9|, we obtain the 
general coupled electron-ion BOA problem 


N -r-,1 N t-7? N 


a= 1 
< n 


a= 1 


a= 1 
\ 


M,y 


+x 

V i=l 

= E k i//k- 


^ _ IT + Vie + Vet 


+ V n 


( 12 ) 


However, one normally separates Eq. (IT2l > into an electronic 
problem only in tp and an ionic problem only in x- To 
do so, one first solves the electronic-only BOA frozen ion 
Schrodinger equation, where the ionic coordinates R„ only 
enter as fixed parameters in tp: 

H e <pi = H e p ( ^'’ R2 ’- RN) (risur 2 s 2 ...r„s n ) 

-£/(R 1 ,R2,...R3v)^ Rl ’ R2 ’- Rw) (riSi,r 2 x 2 ...r ; , i „X 


with 


( 13 ) 
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where 

H e - f e + Vie + Vee + V n . (14) 

In this way, one may find the so-called i th potential energy sur¬ 
faces £,(R 1 ,R 2 , ...R.y) (PES). These are representations of the 
electronic energy landscape as a function of the ionic coordi¬ 
nates. 

In the next step, the ionic BOA Schrodinger equation is 
solved by adding the previously neglected kinetic energy of 
the ions to the potential energy surfaces obtained from the 
frozen ion Schrodinger equation 

i,RiS 2 , --RnSn) = EijXij(R\S i,R2Si, ...R/vSaO, 

(15) 

where 

N 1 

^ = Z-2aT V « + £ ' (Ri ’ R2 ’- R ^ (16) 

a= 1 a 

and the ionic excitations j depend on the electronic excita¬ 
tions, i. 

Comparing Eq. (fl2l) with Eq. (fl6l) we realize that the second 
and third terms of Eq. (Il2l > are neglected in the BOA. This is 
because we assume that the kinetic energy of the ions is not 
affecting the electronic part tp, i.e., V a ip ~ 0. 


C. Ehrenfest dynamics 

Within the Ehrenfest dynamics (ED) scheme fH, we solve 
the coupled evolution of the electrons and ions. The electrons 
evolve quantum mechanically, whereas the ions evolve clas¬ 
sically on a mean time-dependent PES < pj{t) weighted by the 
different BOA PES pi in Eq. i fTH i 

n 

<Pi(Rc(t)) = ^ c i(t)pi (17) 

i 

The ions are evolved according to Newton’s equation of 
motion 

F £D (R ff (0) = (18) 

which satisfies the following potential energy derivative con¬ 

dition 

n 

F* fl (R„(r)) = - 2 lo(0l 2 ( R -«) 

«■ (19) 

= - |v„Ef e (R„(0)| ip(t)j 

where Eq. E3 and the Hellmann-Feynman theorem have 
been employed. The Ehrenfest electron-ion scheme consists 
of the time propagation of the coupled Eqs. O and (O. 


D. Model systems: Initial configurations and Hamiltonians 

We model the positively charged one electron H) and neu¬ 
tral two electron H 2 homonuclear diatomic molecules, assum¬ 
ing their motion is confined to one direction. Such a model 
should provide a reasonable description of a molecule excited 
by a laser field, where the electronic and ionic motion are con¬ 
fined to the polarization axis of the laser field [12]. In this case 
the QMI problem described in Sec. Ill Al where both electrons 
and ions are treated quantum mechanically, can be solved ex¬ 
actly. Furthermore, by working in center of mass coordinates, 
the computational effort required to solve Eq. © is signifi¬ 
cantly reduced. 

However, the singularity in the bare Coulomb interaction of 
Eqs. ©. ©, and <[6]) in ID makes the direct numerical solution 
of the Schrodinger equation © unfeasible. Instead, one em¬ 
ploys the so-called “soft Coulomb interaction” 02,EH]. For 
two particles i and j with charges Q, and Qj , the soft Coulomb 
interaction V m , has the general form 

QiQi 

V int (s ) - , , (20) 

Vs 2 + A 2 

where s is the separation between the two charges and A is 
the soft Coulomb parameter [18;]. Typically, A = ao, although 
other values can also be used |JU)]. 

In essence, the soft Coulomb interaction amounts to a dis¬ 
placement of the trajectories of the two particles in an orthog¬ 
onal direction. So for a hydrogen atom, a soft Coulomb inter¬ 
action of 

V inl (s) =- p- , (21) 

V^ 2+fl o 

is equivalent to having a bare Coulomb interaction with the 
electron and proton trajectories required to be parallel, with 
a minimum separation of ao. This is a quite reasonable as¬ 
sumption, as the most probable electron-proton separation in 
a hydrogen atom is the Bohr radius ao- Soft Coulomb param¬ 
eters correspond to the separation between the ID trajecto¬ 
ries that each electron and ion will move along in 3D with a 
bare Coulomb interaction. We may directly map the ID soft 
Coulomb problem to a bare Coulomb problem in 3D where 
the electrons and ions are separated by the soft Coulomb pa¬ 
rameter distances shown in Fig© while their motion is con¬ 
fined in one direction. In this way, one clearly sees that con¬ 
straint rotations of the molecules are possible in 3D, while 
their motion is still confined to ID. As discussed in Sec. II 11 B] 
the molecules will be perturbed by a kick confined in one di¬ 
rection. 

In Fig© we show the various configurations we have em¬ 
ployed to model an H} or H 2 molecule whose electronic and 
ionic motion is confined to one direction. These configura¬ 
tions are specified by the soft Coulomb parameters between 
the ions A//, the electrons A ee and the ions and electrons A/,,. 
One such configuration has been used previously JH to study 
the dynamics of a one-dimensional H 2 model molecule in 
strong laser fields by means of QMI. 
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MV\ + ^MV\ + ~v 2 + 
1 

1 2 

2 Vl 

1 

P.-Xtf + Al y 

J(X 2 - x 2 ) 2 + A 2 e 

1 

1 

yj(X 2 - Xi) 2 + Aj e A 

j{X 2 - Xl ) 2 + A 2 Ie 

1 

1 

-+ - 



tJ(X 2 -Xi) 2 + A 2 V( x 2 - h) 2 + A 2 e 


respectively. Here, M is the ion mass; Vi, V 2 , X\ and X 2 are 
the ionic velocities and positions for both molecules along the 
direction of motion along the direction of motion. 

The first three and four terms of Eqs. (l22l > and (l23l > are the 
kinetic energies of the electrons and ions in the molecules, as 
explained above. The remaining terms correspond to the at¬ 
tractive and repulsive electrostatic potential energy terms be¬ 
tween such electrons and ions. 

The spatial configuration of positively charged or neutral 
homogeneous diatomic molecules in ID does not change if 
the particle positions are translated uniformly. This reduces 
our three- and four-body coordinate problems into two- and 
three-body ones, respectively. 

We rewrite the classical energies in Eqs. d22l) and (1231 ) 
in terms of the center-of-mass transformation |20|] (see Ap¬ 
pendix [A} to obtain the following two-body ( X,£) and three- 
body ( X,x ,£) Hamiltonians 


FIG. 1. (Color online) Schematic representation of the (a-d) Hi and (e-i) 
H 2 geometries for the minimum ionic separation A//, ion-electron separation 
A i e , and electron-electron separation \ ee for each configuration. Protons are 
shown in red and electrons in black. 


1 d 2 2 M + \ d 2 1 

- AM 3 .« ^ (f+ 0 2 + A 2 

1 1 

V(f-« 2+A /c V* 2+A "' 


(24) 


and 


An analysis of the effect of the initial configuration on 
the optical spectra is shown in Sec. [IV] The classical ener¬ 
gies of positively charged and neutral homonuclear diatomic 
molecules whose motion is confined to one direction are given 
by 


1 ,1 , 1 

E = -MV + -MVl + -vf - 
2 1 2 2 2 


^(v, -X,) 2 + A; 


( 22 ) 


^(X 2 -x l ) 2 + A 2 e y j(X 2 -X l ) 2 + A 2 I 


ll(X,x,A) 


Id 2 d 2 1 +M3 2 
MdX 2 ~ dx 2 ~ AM ~df 

1 1 

^(f-f+^ + A’ V(f - f - ^) 2 + A 2 

1 1 

y](l + l+& 2 + A l 7 ( i + 2-^ )2 + A | 

1 1 

V x2 + A iL ^X 2 + A 2 

(25) 


for positively charged and neutral homogeneous diatomic 
molecules, respectively, after removing the center of mass 
term. Here X and x are the ionic and electronic separations 
and £ is the separation between ionic and electronic centers of 
mass along the direction in which their motion is confined. 
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Although the electrons are treated quantum mechanically 
along their direction of motion, the confinement of their mo¬ 
tion and position along one direction is inherently classical. 
For this reason, our treatment herein is essentially semiclassi- 
cal: quantum mechanical along the direction of motion, and 
classical perpendicular to the direction of motion. This has 
important repercussions for the FF configurations shown in 
Fig- Dig) and (h). For these cases, the Hamiltonian of Eq. ( 1251 
is no longer symmetric under ion or electron exchange, since 
A /( , is either or 2« 0 . This reflects the limitations of such 
a semiclassical treatment. So although the Hamiltonian of 
Eq. (l25l > is still symmetric under ion or electron exchange for 
the H 2 configurations shown in Fig.QJe), (f), and (i), the con¬ 
finement of the electron’s position perpendicular to its motion 
may still have an important impact for these configurations. 

Our aim here is to assess the accuracy of the approxima¬ 
tions introduced in Secs. Ill Bl and III Cl in which the ions are 
treated classically. To accomplish this, we vary the ionic mass 
M in our homogeneous diatomic molecules for the many-body 
problem, while fixing the the ionic charge Q — e. We only 
consider Q - e because the repulsion between the ions of 
more massive homonuclear diatomic molecules with a single 
electron would be so large that the molecules would be un¬ 
stable 0. Furthermore, this allows us to directly compare 
absorption spectra between these model systems for a fixed 
interaction potential. 

E. Symmetries of the many-body wavefunction 

Since for our positively charged homogeneous diatomic 
molecules there are one electron and two ions, the antisym¬ 
metry of the many-body wavefunction must be enforced for 
the ions only as 

MX\S 1 , X 2 S 2 , xs) = -i/f(X 2 S 2 ,X 1 Si,xs), (26) 

for the triplet and 

ifj(X x S \,X 2 S 2 ,xs) = i//(X 2 S2,X l Si,xs), (27) 
for the singlet. 

For our neutral homogeneous diatomic molecules there are 
two electrons and two ions. Therefore, the antisymmetry of 
the many-body wavefunction must be enforced both for the 
ions and the electrons as 

if/(XiSi,X 2 S 2 ,x\si,x 2 s 2 ) - -i//(X 2 S2,X\SuX\S\,x 2 s 2 ), 

(28) 

for the ionic triplet, 

t(X l SuX 2 S 2 ,x l sux 2 s 2 ) = i//(X 2 S 2 ,XiS 1,xisi,x 2 s 2 ), (29) 
for the ionic singlet, 

ij/{X\S\,X 2 S 2 ,xiS\,x 2 s 2 ) = -if/(XiS i,X 2 S 2 , x 2 s 2 , xisi), 

(30) 

for the electronic triplet, and 

\]/(X\S\,X 2 S 2 ,xiSi,x 2 s 2 ) - -t]/(X\S\,X 2 S 2 ,X2S 2 ,x\Si), 

(31) 


for the electronic singlet, respectively. 

Therefore, due to the exchange symmetry of the many-body 
wavefunction, in order to have a total antisymmetric many- 
body wavefunction, the spatial part of the ionic and electronic 
wavefunction must be odd for the triplet and even for the sin¬ 
glet under the exchange of two identical particles. Conse¬ 
quently, we will only be concerned with the spatial part of 
the wavefunction, with the spin part already being separated 
off due to the exchange symmetry of the many-body wave- 
function. 


III. METHODOLOGY 

A. Ground state 

The QMI eigenvalues are obtained by inserting Eqs. (l24ll 
and (l25l > into Eq. (0 for the H) and H 2 molecules, respec¬ 
tively. 

To obtain the PES within the BOA and ED we would insert 
Eqs. (l24l > and ( |25| >. neglecting the first term, into Eq. ( IT3| > for 
the Hi and H 2 molecules, respectively. For the BOA and ED 
ground state electron-ion level, we do not compute Eq. (IT5l) . 
Instead, we fit the ground state PES around its minimum en¬ 
ergy at the inter-ionic distance X eq using a harmonic approxi¬ 
mation E gs (X eq ) + \k\(X - X eq ) 2 , where k\ = <jS[h p is the har¬ 
monic constant, ojj is the harmonic oscillator vibrational fre¬ 
quency and fi p is the ionic reduced mass defined in Eq. (IA5I) . 
From u>i, we obtain the ground state electron-ion eigenvalue 
of a harmonic oscillator s^ )A/ED = E gs (X eq ) + jO>i in the BOA 
and ED PES picture. 

The inversion symmetry with respect to the inter-ionic X 
coordinate of the potential in Eqs. (l24l > and < f25l ). leads to a 
doubly-degenerate solution s\ for each state tl/t in Eq. 0, for 
sufficiently bound global ground state potentials. The inver¬ 
sion symmetry with respect to the inter-electronic x coordi¬ 
nate of the potential in Eqs. (l24l) and (l25l > is not related to 
the statistics of the ions, but to the symmetry of the electronic 
molecular orbital. 


B. Time dependent linear response spectra 

To obtain the linear response photoabsorption spectra we 
apply an initial impulsive perturbation, or “kick” 112111 

7C(Ht) = e iK ( x i +x 2-*) 

(32) 

<K( H 2 ) = e lK<x ' +x ^ x '-^\ ' 

to the ground state wavefunctions i// gs of our H( and H 2 
molecules, respectively, for the BOA and QMI approaches. 
K is a measure of the strength of the kick. We employ a con¬ 
verged kick strength of K - 0.001, for which the linear re¬ 
sponse spectra does not change if it is decreased further. Us¬ 
ing the center of mass coordinates defined in Appendix[A] the 
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terms in Eq. (l32l > become 


for the Hj molecule and 


7C(H 2 + ) = 

7C(H 2 ) = e _iJE2f , 


(33) 


l<A(0> ~ e -iK^e Wl ' (t fo \2f\ %) |^*>. (40) 

k 


where Acivn is the global center of mass coordinate, and £ 
is the separation between the ionic and electronic centers of 
mass along the direction in which their motion is confined. 
The perturbative kick 70 will only induce polarization on the 
coordinates £ defined for H( and Hi in Eqs. (IA3b and (IA8b for 
the BOA and QMI methods. 

In linear response, we expand Eq. ( |33| > in terms of K, ne¬ 
glecting higher order terms 


*(HJ) *1 

7f(H 2 ) » 1 - iK2£. 


(34) 


for the Hi molecule. From Eqs. ( [39b and ( l40b . we see that only 
the ilj gs to odd k dipole moment matrix elements are non-zero 
by symmetry, i.e parity, since f is an odd operator. 

Eq. ( |38| > can be written as 

d(t) ~ ~ 2K 2m + 1 Zi sin 0Jkt \{^ k l^l fa)| > ( 41 ) 

k 

for the Hi molecules using Eq. d39l ) and 

d(t) ~ -4 K Y sin oj k t\(ifr k ||| %)| , (42) 

k 


For the ED approach one should follow the same procedure 
starting from Eq. (l32l ). but substituting the electronic coordi¬ 
nates x for ~ 2 m+ 2 ] £ ~ when M » 1 for H.J and x 2 + x\ 
for —2ij for H 2 . During the time propagation the ions are not 
kicked, but evolve as parameters according to Eq. ( fl9l >. In this 
case, the electron is kicked relative to the center of mass of 
the ions for the Hi molecule, and the two electrons are kicked 
relative to their distance to the ions for the H 2 molecule. How¬ 
ever, the linear response absorption spectra does not depend 
on uniform translations of the ions and electrons. 

The enforced time-reversal symmetry evolution operator 
I 22 ]] we apply to propagate our equations after this external 
perturbation has been applied is given by 

U(t + At, t ) = e -'t«( f+A ') x (35) 

where the Hamiltonian flit + At) is calculated from 


for the H 2 molecules using Eq. (l40l) . where d(t) depends lin¬ 
early on K and u> k = e k - e gs . 

However, in the ED approach the ionic coordinates are up¬ 
dated at each time step. This makes £ and the Hamiltonian 
time dependent. For this reason, the dipole moment from the 
ED approach does not necessarily have the form of Eqs. 6B 
and (l42l >. As we will show in Sec. IIV Al the time-dependent 
effects of zero-point motion within ED, which are incorpo¬ 
rated into the coordinate £(t) = x - Acvi . have an important 
impact on the spectra. 

The optical photoabsorption cross section spectra cr a bs is 
obtained by performing a discrete Fourier transform of d{t) 
(23]]. More precisely. 


Cabs = 47rau>lm 


T 

^YAte^'f^ym-dm 

t =0 


(43) 


>Kt + AO - (36) where 


and the kicked initial state we propagate is 

<A(A0 = e~ iAtH °'K£/g S , (37) 

where i// gs is the ground state eigenstate of the time indepen¬ 
dent Hamiltonian Hq of Eqs. (l24l > and (l25l > for H) and Hi , 
respectively. 

The time dependent Hamiltonian H(t) is then obtained by 
time propagation at each time step self consistently accord¬ 
ing to Eq. (l36l >. starting from the kicked initial state given in 
Eq. ( l37l >. The expectation value of the dipole moment d at 
time t is 

d(t) = <<K0l!hK0>. (38) 

If we assume that the Hamiltonian does not evolve in time 
and we insert Eq. ( l34t into Eq. ( |371 ). using the completeness 
relation Yj k = 1 an d Eq. ([8} we get 


fix) - e ~ 25 *\ (44) 

is a Gaussian damping applied to improve the resolution of the 
photoabsorption peaks, a> is the frequency of the oscillations 
of d(t), a is the fine structure constant, T = 1000 is the total 
propagation time, and At = 0.01 is the time step. 


C. Computational details 

All numerical calculations have been performed using the 
real space electronic structure code Octopus 0. We dis¬ 
cretize the configuration space of the Hi and H 2 molecules, 
using a finite set of values (i.e. a so-called grid) for the co¬ 
ordinates X, x and £ in the box intervals X e [-Lx,Lx), 
x e \—L x , L x ] and £ e [~L C , Lf. These are discretized as 


I <ff(t)) 


'lfc> - iK V 


-is k t 


k 


fk 


2M + 2. 


2M+ 1 




gs 


I >/fk) 
(39) 


Xj = -L x + iAX for i = 0,1,2 ...N x , 
Xj = —L x + jAx for j = 0,1,2 ...N x , 
& = -L f + kAf for k = 0,1,2...% 


(45) 
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using Nx, N x and /V f equally spaced points, respectively. The 
spacing between two adjacent points in the X, x and £ direc¬ 
tions are AX = Ax - A% — -ff-. Convergence is 
achieved when a decrease in AX, Ax, A£ and an increase in 
Lx, L x , L<: does not change the electron-ion static and time 
propagation linear response spectra. 

For the Hj type molecules, ground state convergence is 
achieved for Lx — L( = 1 Oao, AX = 0.05an and A£ = O.lao- 
To obtain the PES we have used Le = lOOao and A£ = O.lao- 
Generally, the convergence of the QMI optical spectra re¬ 
quires L x = 30ao, L{ = 80ao, AX = 0.01 ao and = 0.5ao- 
However, for the ionic mass M of q case, convergence re¬ 
quired Lx - 100ao, L(: - 80«o, AX = 0.03«o and A£ = 0.5flo- 
Finally, for the ED and BOA optical spectra we have used 

Lj: = 500«o and A£ = O.lao- 

For the H 2 type molecules, ground state convergence is 
achieved for L x — L^ = L x = 10ao, AX = 0.07 ao, A£ — 0.2«o 
and Ajt = 0.5ao- To obtain the PES we have used L; - L x — 
40t/() and A£ - Ax - 0.2a 0 . The convergence of the QMI 
optical spectra requires L x = 1 Oa 0 , L x = 80a ( ), L. : = 35«o> 
AX = 0.07«o, Ajt = 0.5flo and Af = 0.6ao. Finally, for the 
ED and BOA optical spectra we have used Q- = L, = 200c/o, 
At; — Ax — 0.5ao. 

Within the BOA and ED, the X coordinate does not need 
to be discretized quantum mechanically. It is either fixed as 
a parameter in BOA or it changes according to the dynamic 
equations in ED. As a consequence, the two and three variable 
bare Coulomb QMI problems confined to ID trajectories for 
the Hj and H 2 molecules respectively, become one and two 
variable BOA and ED problems. These are easier to compute 
numerically, thus providing a more attractive alternative. 


IV. RESULTS AND DISCUSSION 
A. Hf and H> results 

In Fig. [2] we show how the Hj and H 2 BOA ground state 

PES change as a function of the ionic separation yjx 2 + A~ ; 
for each configuration shown in Fig. [I] The PES fitted min¬ 
imum energies at X eq , Eo(X eq ) and positions yjx 2 q + Aj ; are 
shown in Table [j] for the Hj and H2 molecules with the con¬ 
figurations shown in Fig. |T] 

The ground state PESs (dotted black lines in Fig. [2] and 
taken from Ref. |25|) have been obtained by solving the sta¬ 
tionary Schrodinger equation in 3D using basis sets within the 
BOA. Here, the electronic and ionic positions were allowed to 
vary in all spatial directions. The overall shape of these 3D 
PES is reproduced qualitatively by configurations (b) and (c) 
for Hj and (g) for H 2 from Fig. |T] 

The experimental bond lengths of H j and H 2 are 2ao [5?>| 
and V2flo lf25ll . respectively. The equilibrium distance is best 
reproduced by configuration (d) for Hj and (g) for H 2 from 

Fig □ 

The Hamiltonian for configuration (g) for H 2 in Fig.|T|(b) is 
not invariant under electron exchange. Yet, we still consider 



FIG. 2. (Color online) BOA ground state PESs relative to E a (X eq ) 
in eV, versus yjx 2 + Aj ; in A for (a) Ht and (b) H 2 molecules for the 
An, A ee and A/ e configurations shown as insets. The 3D ground state 
PESs (dotted lines) have been taken from Ref.l25l. 

TABLE I. H, and IF ground state PES fitted ground state ener¬ 


gies Eo(X eq ) and positions yX^ + obtained from a harmonic fit 
around X eq for the configurations shown in Fig. |T] 


Species 

An 

Ale 

^ee 

E 0 (X eq ) 

V^ +A « 


(«o) 

(«o) 

(flo) 

(eV) 

(A) 


1 

0.5 

— 

-45.757 

0.5627 

Hj 

0.5 

1 

— 

-21.431 

1.3510 


1 

1 

— 

-21.969 

1.2697 


2 

1 

— 

26.759 

1.0584 


1 

3 

1 . 2 

3 ’ 3 

1 

-60.022 

0.7146 

Hi 

V3 

1 

1 

-45.193 

0.9166 


1 

1 

V3 

-44.790 

0.9004 


V2 

1 

V2 

-45.856 

0.8241 


this configuration, as the H 2 configurations which are invari¬ 
ant under electron exchange (Fig. |TJe,f,i)), yield PES that dif¬ 
fer qualitatively from the 3D PES, as shown in Fig. 0(b). 

The ions sometimes undergo a strong inter-ionic repulsion 
for larger and small X, depending on the initial configuration. 
For the strongly repulsive configurations for small X, the ions 
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FIG. 3. (Color online) Optical spectra for the (a) HJ and (b) Ho 
molecules obtained by classically fixing the ions to their equilibrium 
positions (BOA) and evolving the ions (ED) for masses Mh and M p x 
10 4 with minimum ionic separations (a) A;/ = ao (b) An = A ee = 
V2ao and electron-ion separations (a,b) A= ag shown as insets. 
Evolution of the difference in dipole moment Ad between ED for 
Mh and BOA in milliDebye is shown as an inset of (a). 


are farther apart because the repulsion between the ions is 
stronger than the attraction between the ions and electrons. 
For the strongly attractive configurations for larger X, the ions 
are closer together because the repulsion between the ions is 
weaker than the attraction between the ions and electrons. For 
the latter configurations, more energy is required to dissociate 
the molecule. 

For H P , when A/,, = a o, the potential becomes less repulsive 
for small X as A// increases. However, when Aj e = Q.5ao, the 
potential becomes strongly attractive for larger X. 

For H 2 , the potential becomes strongly repulsive for small 
X for the linear configuration (g). When A// = ao and A j e 4- ao, 
configuration (d) in Fig. Q] the ground state PES is unbound. 
As the electrons are necessarily very close to each other (A,,,, = 
|ao) when the molecule is bound, their repulsion forces the 
dissociation of the H 2 molecule into two isolated stabler H 
atoms. For this reason we will disregard this configuration 
from hereon. 

In Fig. [3] we compare the H-J and H 2 optical spectra we ob¬ 
tain by classically fixing and letting the ions evolve according 
to ED in time from X eq . Essentially, including the classical 
movement of the ions hardly changes the spectra. However, 
new peaks appear before the first electronic excitation for both 



10 20 
tins (eV) 


FIG. 4. (Color online) Absorption spectra obtained from a clas¬ 
sical BOA (dashed lines) or quantum mechanical QMI (solid lines) 
treatment of the ions of an (a-d) HJ molecule with configurations 
(a) A„ = \a 0 ',Ai e = a 0 (green), (b) A„ = 2a 0 ;A /e = a 0 (blue), (c) 
An = ao; A h = ao (red), and (d) A /; = a 0 ; A ie = |a 0 (orange) or an (e- 
h) Hi molecule with configurations (e) A u = a 0 ; A /f = ao; A ee = lag 
(violet), (f) A„ = V2a 0 ; A Ie = a 0 ; A ee = V2a 0 (blue), (g) A„ = 
V3ao; A Ie = ao ; A ee = a 0 (red), and (h) A« = |a 0 ; A/ e = |, |ao; A ee = 
a 0 (green) shown as insets. Dotted vertical lines denote the energies 
Si of the unoccupied electronic levels ip j (£) relative to the ground 
state energy eo for each configuration at X eq . 


the and H 2 molecule, at 1 and 12 eV, respectively. For H 4 , 
the new peak corresponds to the frequency of the ionic zero- 
point motion around X eq , which vanishes for large masses 
(M p x 10 4 ) because heavy ions hardly move around X eq . On 
the other hand, H 2 ’s higher energy peak does not vanish for 
large M. Due to its width, the ED and fixed ion spectra do 
not overlap completely. We explain the origin of this peak in 
Sec. HVCl 

The inset of Fig. [3 a) illustrates the time-dependent effects 
of zero-point motion. In the BOA the ionic center of mass 
Xcmi is fixed, so the electron can only oscillate about it. ED 
(and QMI), however, allow ionic motion, so long as the global 
center of mass Acm, is conserved. The increasing difference 
Ad between the ED and BOA dipole moments thus demon- 























































































































strates the ions move, e.g., Ac/(24fs) ~ 0.7mD. 

In Fig. |4j we show how a quantum mechanical treatment 
of the ions (QMI) affects the optical absorption spectra for 
HI and FF molecules in the configurations of A//, A,,,, and A/ f 
shown in Fig. [T] We see that new features emerge in the spec¬ 
tra when the ions are treated quantum mechanically instead 
of classically. The peaks are broadened, become asymmetric, 
and their amplitudes and energies change as a function of the 
initial configuration and charge of the molecule. In particular, 
comparing the BOA and QMI spectra shown in Fig. [4] we find 
that each peak splits into a lower and higher energy contribu¬ 
tion. Depending on the energy shift and amplitude of each 
contribution, these can appear as separate peaks or shoulders 
in the spectra. The shoulders are giving rise to an asymme¬ 
try that can be seen for almost every peak. These quantum 
features are not as strong for the neutral FF homonuclear di¬ 
atomic molecule, regardless of the initial configuration. With 
a classical description of the ions, we do not obtain these 
quantum mechanical features in the optical spectra. 

Generally, we find treating the ions quantum mechanically 
substantially affects both the peak positions and widths in the 
absorption spectra for most of the configurations considered. 
For inter-ionic potentials which are less repulsive (Fig. [4}b) 
and (e)), the line shape of the QMI peaks is narrowed, and 
approaches the fixed-ion at X eq limit. For potentials which are 
attractive for larger X (Fig. [Jjd) and (f)), all the QMI peaks 
are blue shifted with respect to the fixed-ion at X eq spectra. In 
this case, the peak excitation energies are larger because more 
energy is required to excite these transitions. 


B. Mass dependency 

To provide a quantitative analysis of the differences be¬ 
tween a classical (BOA/ED) or quantum (QMI) treatment of 
the ions, we will compare the total ground state energies and 
the peak positions and widths in the absorption spectra as we 
vary the ionic mass over seven orders of magnitude. 

The accuracy of the static BOA and ED calculations can be 
understood from a perturbation theory argument in terms of 
the small parameter k — (m e /M) /t , defined as the ratio between 
the ionic and electronic displacement (25$], where X — X eq + 
k(. We illustrate this in detail in Appendix [B] 

To test the accuracy of the BOA and ED approximations, 
we first compare the BOA and ED ground state electron-ion 
eigenvalues to those obtained from QMI. We expect that the 
BOA and ED should be accurate around the minimum of the 
ground state PES, as the exact eigenvalues of the electron-ion 
problem can be interpreted in terms of the ionic vibrational 
levels for the electronic ground state PES. As discussed in 
Sec. IIII Al the ionic contribution comes from the ground state 
level of a quantum harmonic oscillator where the mass is in¬ 
cluded via a?/. 

The ground state electron-ion eigenvalue for H) and PE 
molecules whose motion is confined in one direction is given 
by 

(46) 



FIG. 5. (Color online) Difference in ground state total energy be¬ 
tween QMI and BOA/ED approaches As gs in eV versus the electron- 
ion mass ratio m e /M for HJ (A;; = A Ie = ao; O) and FE (A// = A fc = 
V2ao and A = «o; ■). Solid lines are a power law fit a {m e /M) b . 


TABLE II. M) ground state eigenvalues obtained from diagonaliza- 
tion of the QMI approach s g S MI , and ground state harmonic BOA and 
ED £^ j ° 4/£d eigenvalues obtained from a harmonic fit around the min¬ 
imum of the ground state PES. We show these results for different 
ionic masses M and A n = A Ie = a 0 . 


M 

4 M (eV) 

s * 0A/ED (eV) 

m e 

M 

P 

-21.703395 

-21.697297 

0.004836 

H 

-21.970225 

-21.969323 

0.000545 

D 

-22.009642 

-22.009272 

0.000272 

T 

-22.027041 

-22.026787 

0.000182 

Li 

-22.053541 

-22.053420 

0.000079 

Na 

-22.076626 

-22.076580 

0.000023 

K 

-22.083189 

-22.083160 

0.000014 


where e (0) and s i2> correspond to the electronic and ionic mo¬ 
tion eigenvalues, and s (1) and e® are equal to zero by sym¬ 
metry (see Appendix [B]). The first order correction to the 
ground state energy for the full electron-ion problem using 
the BOA/ED is the term of fourth order in k. 

To check the dependence of the static ground state eigen¬ 
value accuracy of the BOA and ED approaches, on the 
electron-ion mass ratio, we use the following power law re¬ 
lation 

(47) 

Note that most of the molecules used in this analysis are fic¬ 
titious because we do not change the charge of the ions as 
explained in Sec. MDI except for H, D, and T, as these have a 
positive electric charge of e. 

Fitting the ground state error from Eq. (l47l > to the data in 
TablesHIlandllin we obtain a power law of/? ~ 0.92(2) and b ~ 
1.05(2) for the BOA/ED approaches, as shown in Fig. [5] This 
means the BOA and ED energy expression gives the correct 
total ground state energy of the full electron-ion problem up 


4 0A W°) + e ( 2 V) + 0(K 4 ), 
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TABLE III. M 2 ground state eigenvalues obtained by electron-ion 
QMI diagonalization £g S MI , and ground state harmonic BOA and ED 
gBOA/ED vibration levels obtained from a harmonic fit around the min¬ 
imum of the ground state PES. We show these results for different 
electron-ion mass ratios ^ and A ;; = V2flo and A= Go¬ 


M 

4 MI (eV) 

(eV) 

m e 

M 

M- 

-45.621166 

-45.440564 

0.004836 

H 

-45.739425 

-45.722955 

0.000545 

D 

-45.772460 

-45.763130 

0.000272 

T 

-45.786936 

-45.779126 

0.000182 

Li 

-45.810176 

-45.807442 

0.000079 

Na 

-45.831319 

-45.830529 

0.000023 

K 

-45.837552 

-45.837070 

0.000014 


to fourth order in n. 

In Figs. [ 6 ] and [7] we show how the absorption spectra de¬ 
pends on the mass for the H) and FL configurations for which 
the overall PES shape is closest to that from the 3D treatment 
in Ref. [25l Specifically, we analyze the Hi configurations 
shown in Figs. |TJb) and (c), and the FL configuration shown 
in Fig. [Hg). 

In Fig. Ob) we see that in the large mass limit ( M ~ 10 4 x 
M p ), the QMI spectra exhibits even to odd transitions which 
are allowed by the symmetry of the electronic wavefunctions 
shown as insets. For every allowed transition in Figs. [ 6 ] 
and [7] we have a red shifted and a blue shifted contribution. 

The position of the first, second and fourth peaks in 
Fig- Eta) and (b) (to 1 , u> 2 , and 0 J 4 ) are red-shifted and the third 
and fifth peaks (to 3 and 005 ) are blue-shifted with respect to the 
fixed-ion at X eq spectra. As the mass increases, all peaks tend 
towards the fixed-ion at X eq limit. In Fig. [7] all the peaks are 
red-shifted, although the second peak is also a classical peak 
as shown in Fig. [4jh), which disappears for smaller masses. 

In Fig. [ 8 ] we show the symmetry of the occupied electronic 
wavefunction tp 0 (i ;, x) and the first ten unoccupied electronic 
wavefunctions 1 p,(^, x) for H 2 . Only transitions to unoccupied 
electronic wavefunctions that are even functions of x and odd 
functions of 'c should contribute to the absorption spectra by 
symmetry, i.e. (< p ,(^, jc)|| |^>o(£, *)) > 0. However, Fig. Hfh) 
shows there is an absorption peak in the BOA spectra for the 

x 

—» <p 3 transition, despite ip 3 eq (^,x) being an odd function 
of x, as shown in Fig. [ 8 ] This is because the Hamiltonian for 
the configuration A// = Ifl 0 , A« = ao , and A j e = |«o is n °t 

invariant under electron exchange. 

From the ED spectra in Fig. Ob), we also have an additional 
peak at a lower energy. When the ions are fixed, this peak 
is less intense than when they are allowed to evolve. In this 
case, the peaks’ energy is given by the first excited transition 
((fio —> <f\ ) as seen from the energy of the vertical dotted frozen 
ion lines in Fig-IHf). However, the unoccupied wavefunction 
ip \(£, x) is even with respect to 'c and odd with respect to x, 
as shown in Fig. [ 8 ] This suggests such a transition should 
initially be parity forbidden. 

To calculate the dipole moment for H 2 and different M, 
we only kick our molecules along <f, as shown in Sec. IIIIBl 



0 5 10 15 20 25 30 35 40 


ha (eV) 


FIG. 6. (Color online) QMI electron-ion absorption spectra for a 
positively charged homonuclear diatomic molecule with ionic mass 
M of e, p., H, D, T, Li, Na, K, or 10 4 p in the configuration (a) A/; = 
^ao',Ai e = ao or (b) A /; = A Ie = ao. Dotted vertical lines denote 

X 

the energies e, of the unoccupied electronic levels ip j (Q (shown in 
blue as insets) relative to the energy £ 0 of the ground state electronic 
level <fi 0 rq (^) (shown in green as insets) for each configuration. Note 
that the spectra have been scaled with decreasing mass for clarity. 
Portions of (b) have been adapted from Ref.l30l. 



ha (eV) 


FIG. 7. (Color online) QMI electron-ion absorption spectra for 
a neutral homonuclear diatomic molecule with ionic mass M of H, 
D, T, or Li in the configuration A n = ~a 0 , A ee = ao, and A le = 
|flo- Dotted vertical lines denote the energies £,- of the unoccupied 
electronic levels ip*" 11 (i;,x) (shown as insets) relative to the ground 
state energy £ 0 for each configuration. Note that the spectra have 
been scaled with decreasing mass for clarity. 
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x 

FIG. 8. (Color online) Electronic wave functions <p j eq (g,x) for i = 
0,..., 10 of an H 2 molecule in the configuration A« = ja 0 ; A ee = a Q . 


When calculating the spectra in the BOA, the ionic coordinate 
is frozen at X eq and cannot evolve in time. The electronic co¬ 
ordinate x forms part of the integral, but can evolve in time. 
When we apply a kick along the distribution of the charge 
in the molecule will change with time. The electrons and ions 
will feel the charge distribution of the other particles. Thus, 
the electronic coordinate x can evolve in time, although this 
effect is not taken into account when the dipole moment is 
calculated. Essentially, the (£, x) basis is rotated by the kick 
to a (£', x') basis. 

As £ is time dependent within ED, the tp\(£', x') rotated ba¬ 
sis has a mixture of even and odd components in both x and £, 
removing the parity constraint on the tpo —» q>\ transition. 

As shown in Fig. [Tjb), this extra parity forbidden peak is 
not as intense as the other peaks which are allowed by sym¬ 
metry. This peak is rather weak because the electrons and ions 
are close to each other, but not on the same plane, for the con¬ 
figuration shown in Fig. [lji). Additionally, as the ipo —> ip 1 
transition is initially forbidden by symmetry for both x and <f, 
the rotated contribution with even and odd symmetry in x and 
£ is small. 

Overall, heavier ions have narrower peaks as we approach 
the classical limit. Figure [9] presents this effect in the time 
domain. Energy transfer between the excited electrons and 
the ionic system is already clearly seen after a few fs. Even 
in the large-mass limit (M p x 10 4 ) energy transfer is clearly 
evident. Oscillatory behavior, including beat frequencies, is 
still present 24 fs after the initial kick. 

When the ions evolve quantum mechanically, the electrons 
can transfer part of their dipole moment to the ions. The am¬ 
plitude of the dipole moment thus decreases at different rates 
depending on the ionic mass. This process will take longer as 
the mass of the ions increases and it becomes more difficult 
to displace the ions. For very large ion masses, the interac¬ 
tion with the electronic motion becomes nearly elastic. This 
allows the electrons to oscillate back and forth without the in¬ 
fluence of any external ionic displacements. Since the widths 
of the absorption peaks are proportional to the energy transfer 
from the electrons to the ions, we expect the widths to scale 
as the electron-ion mass ratio to the one fourth, as discussed 
in Section lHBl and AnnendixlBl 



FIG. 9. (Color online) Evolution of the QMI electron-ion dipole mo¬ 
ment d in milliDebye with time t after the initial “kick” in fs for ionic 
masses M of p, H, D, Ti, Li, Na, K, or (upper inset) 10 4 p of a pos¬ 
itively charged homonuclear diatomic molecule in the configuration 
An = A /e = ao . (lower inset) After 20 fs, the amplitude has decreased 
by a factor of one million. 

C. Spectral lineshape 

To quantify the width and energy of the peaks in the spectra, 
we have employed both Gaussian 

X (48) 

1=1 

and Lorentzian 

V/, -(49, 

functions. Here /, is the intensity, w, the position, cr, the stan¬ 
dard deviation, and T, the full width at half maximum of the 
first three peaks of the QMI spectra. 

From Fig. [TO] in which we show the QMI spectra for fQ, 
we clearly see that the tails of the peaks of the QMI spectra 
are Gaussian. Moreover, the three peaks can only be fitted si¬ 
multaneously with Gaussian functions, as the Lorentzian fit to 
the first peak decays so slowly that the second and third peaks 
are completely obscured. Furthermore, the ionic wave packet 
on the ground state PES is a solution of a harmonic eigen¬ 
value problem and thus should have a Gaussian line shape. 
This means the spectral line shape arises from the shape of 
the PES, rather than the coupling between ionic vibrations of 
the molecule. 

Note that the width of the fixed-ion at X eq spectra in Fig. [3] 
is due to the artificial damping introduced in the spectra. The 
electronic transitions should be delta-like functions, but are 
convoluted with a Gaussian function to plot the spectra (see 
Eqs. ( l43l ) and (l44ll ). However, the widths in the QMI spec¬ 
tra are physical, and the Gaussian line shape is due to the 
electron-electron coupling via the ionic displacements. 
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FIG. 10. (Color online) QMI electron-ion absorption spectra for an 
HJ molecule in the configuration A ;/ = A /c = a 0 obtained with (red 
solid lines) and without (green dashed lines) imposing symmetry in 
X on the ionic wave functions. Lorentzian (blue dotted lines) and 
Gaussian (black solid line) fits to the first three peaks of the spectra 
are also provided. 


To ensure that the optical spectra we obtain is only affected 
by the external perturbation K, we have used a symmeterized 
initial wavefunction 




V2 


(50) 


Thus, we always excite from a ground state which is sym¬ 
metric in the ionic coordinate. In Fig. [TO] we show that sym¬ 
metrizing the wavefunction does not change the calculated op¬ 
tical absorption spectrum. This means that we already obtain 
a nearly symmetric ground state starting configuration from 
the stationary Schrodinger equation ©. However, the data 
shown in Fig.[ 6 jb) has been symmetrized for every M. 


D. Model 

To explain why a quantum treatment of the ions has such a 
strong effect on the absorption spectra for H 5, we propose a 
simple two level model. Using this model, we will show how 
the observed QMI spectral peaks and widths can be extracted 
from the electronic BOA eigenenergies s, at equilibrium X eq 
of the ground state through the electron-ion mass ratio m e IM. 

When an external kick is applied, a charge separation is in¬ 
duced in the molecule which will oscillate back and forth with 
time. As discussed in Sec. IHIBI the applied kick is simply a 
transformation of the ground state wavefunctions, through the 
application of a phase factor e lKt , to eigenstates of the sys¬ 
tem with momentum K along the direction of motion. This 
leads to a transition dipole moment between an initial and a 
final electronic state. 

The electronic wavefunctions with even indices (p 2 i are even 
functions of while the electronic wavefunctions with odd 
indices ifi 2 i+\ are odd functions of £. This parity of the elec¬ 
tronic wavefunctions means that the transition dipole moment 


is zero for transitions from the ground state to even unoccu¬ 
pied states, i.e., (^ 2 ;+ 2 lfl</ J o) = 0. Essentially, optical transi¬ 
tions (po —> ip 2 i +2 are forbidden so long as <f 2 i +2 is an even 
function of <f. 

This is the case when the ions are treated classically. In 
fact, Figs.[4ja-d) clearly show that the peaks in the absorption 
spectra obtained from a classical BOA treatment are always 
aligned with the energies of odd-parity unoccupied electronic 
levels S 2 i+\, i.e., the Franck-Condon transitions ifo —> ip 2 i +1 ■ 

When the ions are treated quantum mechanically, every al¬ 
lowed transition is split into red and blue shifted contributions, 
with the shifts increasing as the mass decreases. The level 
splitting we observe in Fig. [ 6 ]is reminiscent of level hybridiza¬ 
tion. 

This motivates us to employ a simple two-level model ll27l 
l28tl to describe the energies and widths of the QMI peaks. 

To do so, for each odd-parity unoccupied electronic level 
ip 2 i+i at fc' 2 ,+ 1 , we artificially introduce a level at ejj+i to which 
it couples. 

As mentioned in Section IIIB I and Appendix [B] the ratio 
between the vibrational and electronic energies, E v n,/E e i ec , 
scales as the square of the ratio between the ionic and elec¬ 
tronic displacement (6 /ao ) 2 (29]. This means the ionic dis¬ 
placement scales as the electron-ion mass ratio to the one 
fourth 6 ~ (m e /M)' A . Since the probability of coupling is 
directly related to the quantum ionic displacement, we ex¬ 
pect the coupling between the energy levels to scale as 6 w 
(m e /M) v \ Further, the width of the peaks in the absorption 
spectra should also be related to the ionic displacement. We 
thus assume the coupling between the energy levels is propor¬ 
tional to the ionic displacement 6 ~ M i.e, aM 7“ where a 
is the constant of proportionality. 


The resulting two-level Hamiltonian 


S 2/+1 

a 

"W 


solution for the S 2 ,+i allowed peak’s energy of 


a 

aF 

S 2/+1 


has a 


OJ2i,2i+l 


£2i+l + 62/+1 
2 


l 


£2i+l - 62;+1 
2 



(51) 


This two-level model yields two peaks that are lower and 
higher in energy, through level repulsion. As the mass de¬ 
creases, ionic displacements become larger. This leads to a 
greater coupling. The coupling between the energy levels will 
be larger when the mass decreases and the energy separation 
between the two coupling energy levels will also increase. 

In Figs.QT]and[T2]we use the two-level model to fit the cal¬ 
culated QMI peaks for various ionic masses M. Specifically, 
we employ Eqs. (fsTb to fit oq, u> 2 , and aq for H^ and H 2 . In 
general, the calculated peak positions are within 0.1 eV of the 
two-level model fit, which is also the expected accuracy of 
such calculations. In each case, we find the coupling between 
the transitions has a constant of proportionality of a « 11 eV. 
Furthermore, the artificial level S 3 « 23.5 eV for both config¬ 
urations of HJ. 

For H*, we find the first peak depends on M 1/2 , since 
a 2 /M V: <s H\ - si * 9.7 eV. In other words, 

a 2 /M' /2 
«i ~ ei - --- 

si - Si 


( 52 ) 
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FIG. 11. (Color online) Two-level model fits to the first four peaks in the QMI electron-ion absorption spectra <u, for a positively charged 
diatomic molecule with ionic mass M in the configuration A Ie = a 0 and A u = \a 0 (•) or A n = a 0 (O). Level coupling has the form a(m e /My'\ 
and the decoupled levels are obtained from the ground state electronic eigenenergies £,. Gray regions denote a ±0.1 eV estimated accuracy. 




FIG. 12. (Color online) Two-level model fits to the first three peaks 
in the QMI electron-ion absorption spectra cu, for a neutral diatomic 
molecule with ionic mass M in the configuration A/ £ = |; |fl 0 . A;/ = 
|a 0 , and A ee = a 0 (■). Level coupling has the form a(m e /M) v \ and 
the decoupled levels are obtained from the ground state electronic 
eigenenergies £,. Gray regions denote a ±0.1 eV estimated accuracy. 


where ~s\ ~ 20 eV. 

For the second and third peaks, so long as £3 - £3 <s: 
2 a/M'k, we may further approximate the peaks by 


^ 2,3 


S3 +«3 

2 


a (£3 - £ 3 ) 2 M i/j 

W' + 8cr 


£3 + £3 a 

2 * W'’ 


(53) 


As we see in Fig. [IT] this is indeed the case for the second and 
third peaks in the absorption spectra, 102 and M 3 , of H) in the 
configuration A j e = ciq and A// = ~ao, as £3 - £3 « 0.4 eV. 

Essentially, all the peak positions in the QMI spectra are fit 
using only two parameters, the artificial level’s energy £ 21 + 1 , 
and the coupling constant a. However, we are not always able 
to decouple the two peaks’ red and blue shifted contributions 


FIG. 13. (Color online) Full width at half maximum (FWHM) of 
the Gaussian fits (FWHM = 2 V5cr,) to the first and third peaks of 
the absorption spectra for a positively charged homonuclear diatomic 
molecule with ionic mass M of p, H, D, T, Li, Na, or K versus the 
fourth root of the electron-ion mass ratio ( m e /M) v ‘ in the configura¬ 
tion A Ie = a 0 and A n = (•) or Au = a 0 (O). Black lines are 

linear fits to each peak for both configurations. Gray regions denote 
a ±0.1 eV estimated accuracy. 


because of their overlap due to their finite width. This is par¬ 
ticularly true for Hi. As a result, we have fewer data points 
for the H 2 peaks, as shown in Fig. [T21 reducing the reliability 
of the fit to Eq. |5T] 

In Fig. [13] we also show that the width of the first and 
third peaks for the configurations shown in Fig. □ (b) and 
(c) scale as the electron-ion mass ratio to the one fourth, 
i.e. FWHM ~ (m e /M)' / \ as expected from our model. The 
FWHM for the first peak has a larger constant of proportion¬ 
ality than the third peak, but the widths for both configurations 
may be fit simultaneously. Altogether, this demonstrates the 
predictive power of the simple two-level model for describing 
the QMI spectra as a function of the ionic mass. 
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V. CONCLUSIONS 

We have shown that additional features may appear in the 
linear response spectra of charged and neutral H 2 homonu- 
clear diatomic molecules when the ionic motion is described 
quantum mechanically. Such features are strongly dependent 
on the molecules’ configuration, i.e., the shape of the PES. 
The widely used classical ionic motion BOA and ED ap¬ 
proaches fail to describe such features. We also demonstrate 
that these features may be understood using a predictive two- 
level model. These results demonstrate how for light atoms, 
the quantum nature of the ions may play an important role 
when describing absorption processes. 


where Acm 2 is the global center of mass coordinate and £ is 
the distance between the electron x and the ionic center of 
mass Acm, . Here the velocities are the time derivatives of the 
positions. 

By substituting Eqs. ( lAll i. (1 A2 b . and (IA3b into Eq. (1221 ) 
we obtain for the classical energy of a positively charged 
homonuclear diatomic molecule 

E = 1(2 M + l)V* Mi + + X -» e V] 

1 1 1 

Vd + f ) 2 + A i V ( t- 0 2 + a * 

(A4) 
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Appendix A: Center of mass transformation 


Ee = 


2 M 


M 


2M + 1 ’ IUp 2 


(A5) 


If we rewrite Eq. <IA4b in terms of the momenta given by 

d 


Px cm 2 = (2 M + 1)V C m 2 = -i 

Px = » P V X =-i-^. 


9Xcm 2 ’ 


(A6) 




Here we present in detail the coordinate transformations ap¬ 
plied in our description of a homonuclear diatomic molecule 
whose electronic and ionic motion has been confined to one 
direction. First, we perform a center of mass transformation 
of the two ionic coordinates X\ and A3 


_ X\ + X 2 _ _ V\ + V 2 

CMl “ 2 ’ CMl ” 2 ’ (Al) 


X = X 2 - X]\ 


Vx = V 2 -Vu 


we obtain the two-body Hamiltonian in Eq. (l24ll . 


2. Neutral homonuclear diatomic molecule 

For a neutral homonuclear diatomic molecule with two 
electrons, we also perform a center of mass transformation 
of the two electronic coordinates xi and x 2 


where Acm, is the center of mass coordinate of the ions and 
X is the distance between the ions. Here the velocities are the 
time derivatives of the positions. 


X\ + X2 Vi + V2 

*CM, = -^-> l ' CM i = -2-’ 

X-X2-X1; V = V2~Vi, 


(A7) 


1. Positively charged homonuclear diatomic molecule 

For a positively charged homonuclear diatomic molecule 
with one electron, the electronic coordinate and velocity are 
simply 

x — x\; v = v'l. (A2) 

Next, we perform a global center of mass transformation of 
the center of ionic mass and electronic coordinates Acm, and 
x, keeping the ionic separation X fixed 

2MX C Mi +*_ _ 2MV C m, +v 

CM2 ” 2M + 1 ’ CM2 “ 2M+ 1 ’ (A3) 

£ = x - Acm, ; Vf = v - Vcm, , 


where xcm is the center of mass coordinate of the electrons 
and x is the distance between the electrons. Here the velocities 
are the time derivatives of the positions. 

We now perform a global center of mass transformation of 
the two ionic and electronic center of mass coordinates coor¬ 
dinates Acmi and xcm! , keeping the ionic and electronic sepa¬ 
rations A and xfixed 

2MAcm, + 2xcm, _ 2MVCM, + 2vcm, 

* CM2 ~ 2M72 ; yCM2 “ 2MT2 ’ 

f = -xtM, - Acm, ; Vf = vcm, - Vcm, , 

(A8) 

where Acm 2 is the global center of mass coordinate and ii is 
the distance between Acm, and xcm, ■ Here the velocities are 
the time derivatives of the positions. 
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By substituting Eqs. (IA 1 1) . (IA7I) . and (IA 8 b into Eq. (f23l) we 
obtain for the classical energy of a neutral homonuclear di¬ 
atomic molecule 


E =1(2 M + 2)v£ M2 + \p p v\ + \lLe P V] + \]ieV 2 


2 A 2 ^ * 2 

1 1 


^(f-f+02 + A?, ^(f - f - 0 2 + A; 


(A9) 


V ( f+ f + ^) 2 + A / f a/ ( ? + 2-^ 2 + A * 


^X 2 + A 2 V a ' 2 + A L ’ 


where fj. p is the reduced mass of the two ions, is the re¬ 
duced mass of the two ions plus two electron system, and p ( , 
is the reduced mass of the two electrons 


M 


2 M 


E p - 2 ; Eep - l+M 


1 

; a = 2 - 


(A10) 


If we rewrite Eq. i l/\9| i in terms of the momenta given by 

d 


dXcM 2 


(All) 


P XcMl = (2 M + 2)V C m 2 = -i 

P x =„ pVx = -i± 

. d 

~ A l epV( ~ 

- _ d 

Px = EeV = -l —, 
ox 


; obtain the three-body Hamiltonian (X,x,£) in Eq. ( l25t . 


Appendix B: Accuracy of the BOA 

Here we provide a detailed analysis of the accuracy of the 
Born-Oppenheimer Approximation (BOA) El] for the case of 
a homonuclear diatomic molecule whose electronic and ionic 
motion is confined to one direction. In general, the ratio of 
vibrational to electronic energies, EA to E e j ec depends on the 
electron-ion mass ratio m e /M as [31] 


^vib 
Pelec 


(Bl) 


k = ( m e /M )'/* 


to third order as follows: 


*(*, ♦*.*.*)»-A*. + 

+ *—V(X,x,£) f 
OA x=X cq 

id 2 \ 7 d 2 

~x~—x + AT* x^V(X,x,£) 


<9<T 2 '2!" dX 2 
1 , <9 3 

+ A w v(Xx ’ e> 


x=x eq 

? + o(A 


x=x„ 


Sorting the Hamiltonian in different powers of k, i.e., 
H(X eq + <, x, f) * H (0) + x// (1) + k 2 H (2) + k 3 H {3 \ 
we obtain: 

1 d 2 


(B2) 


(B3) 


1 /9 2 

H( 0 ) = _ JL £_ 


2 ^c <9x 2 2 fi e p d £ 2 


+ V(X eq ,x,{). 


H a) = —V(X,x,& 
dX h 


i. 


r ) 2 1 /5 2 

^ C2) = -^ + 2 !^ yK ^ 


(B4) 






1 (9 


<T 3 - 


Expanding the time-independent Schrodinger equation © 
in powers of k to the third order, we obtain: 

3 3 

Y J (K n P (n) )[i< n ^ n) ] = ^(/V' i} )[k> (,,} ]. (B5) 

«=0 n=0 

Decomposing Eq. ( IB5 b in terms of k, we find 

O(/c 0 ) :H (0 V (0) > = £ ( °V (0) >, (B6) 

H (0 V (1) ) + H (1 V (0) > = g (0 V (1) ) + £ (I V (0) >, (B7) 


O (x 1 ) 

Olx 2 ) 


H ( 0 V (2) ) + + H (Z W U> ) 


Wi.z/Dn 


>( 2 )|,/,( 0 ) n 


= e (0) |^ 2) > + £ ( 1 V (1) > + £ ( 2 ) I-A (0) >, 


(B 8 ) 


0 (x J ) :H < 0 V (3) > + H (L, \ifr (z> ) + + H K 3 , \iff w ) 


KD..A2) 


(2)|./AD\ 




= e (0)|,/,(3) 


|(/r (3) > + £ (1 V (2) > + £ (2 V (1) > + £ (3 V (0) >-(B 9 ) 

H a)) is the electronic frozen ion Hamiltonian at X eq and e (0) 
is the zeroth-order eigenvalue which corresponds to the elec¬ 
tronic motion. Therefore, we choose the zeroth-order wave- 
function as: 


^Xx eq ,z,0 =x(0<p m (x eq ,a 


(BIO) 


where 6 is the length scale of vibrational motion, and a () is 
the length scale of electronic motion, i.e., the Bohr radius. 
This means the ratio of ionic to electronic motion is of the 
order 6 /ao ~ ( m e /M With this in mind, we may expand the 
Hamiltonian in Eq. (l25l > as a function of the small parameter 


where ip ii]> is the electronic ground state wavefunction of H a)) 
and x (0 is the ionic wavefunction which will be specified 
later. 

Based on Eqs. ( IB7 b . ( IB 1 Ob and the Hellmann-Feynman the¬ 
orem, A 11 vanishes. This is because the first derivative with 
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respect to the eigenvalue e (0) at X eq is zero. More explicitly, 

£ (1 > = <t/r ( 0 ) |// ( 1 V (0) ) 

^ (0) Wifi*> 


= y (0) 




dX 


P°\X) 


x=x,_ 

<x\Z\x) = o. 


(fill) 


From Eq. (IB 8 b we obtain the second-order correction to the 
energy: 


5 2 


X 


dp 


X 


JO) 


5 F W V(X ■ x 'f> 


<P (0) )(x \e\x) (B12) 


Z 

n>0 


PPUv(X,x,p\ 


(0) 


JO) - J°) 


■<x\?\x). 


where the first order correction to the wavefunction is ob¬ 
tained from Eqs. ( 1B7 | i and ( IB 10b : 




n>0 


■z 

n>0 


JO) (0) 


w?) 


PPUv(X, x p)\LJ^) 


- 1 ^). 


jo) _ jo) 

*0 


(B13) 


Here £, ( , 0) and |^!, 0) ) are the n ,h electronic eigenvalue and eigen¬ 
state of the Hamiltonian H [l)] . 

We now choose x(0 (see Eq. (IB 10b ) to be the lowest eigen¬ 
function of the harmonic oscillator problem. We can then ex¬ 
press t: <2> in the form 


s (2) = \ x 


P_ 

'dp 


X^j + (*" k 2 U) = \<>>u (B14) 


where 


k\ =lp 0) 


dX- 


■V(X,x,P 


JO) 


~ 2 Z 

n>0 


P?\±v{x, x ,o\Ul 


(0) 


J 0 ) JO) 


(B15) 


is the harmonic oscillator constant. Second order corrections 
to the energy thus correspond to the ionic vibrations. 


Finally, from Eqs. ( IB9b . (IB 1 0b and (IB 13b we obtain for the 
third-order correction to the energy: 


e (3) =<t// o) [H (i y (2) > + y°w 2 y (1) > + <i/r (0) i // ( 3 y (0) > 


= U 0) 


— V(X,x,p 
dX ^ 


x=x e , 


v (2) ) <x\F\x) ~ [x 
J ) \ < x\C~\x) 


+ [<p 


J 0 ) 


J0) 


1 d 2 
2 ! 

1 <9 3 


dp 


X 


x=x„ 


V (0) )<X\P\X ), 


(B16) 


where the second-order correction to the wavefunction is from 
Eqs. « and (lBl3b : 


kO=z 


/ 2>0 


(* 2 ®| H^) 

JO) _ JO) 




z 

«>0 


(PP\sV-HV\P*) 


(B17) 


J°) _ J0) 

C'n 


l^ 0) >l^o 0) >- 


All the terms from Eq. (IB 16b using Eqs. (IB 13b and (IB 17b 
vanish by parity. This is because they are all proportional to 
<t/r (0) |// t3 y (0) > which is zero by parity since H (3) is odd in ( 
and p 0> is even in p For this reason e (3) = 0, and the error 
in the BOA ground state energy, after including the zero-point 
energy correction, is 0(/e 4 ) ~ m e /M, as shown in Fig.0 

For the neutral homonuclear diatomic molecule we follow 
the same procedure as above using Eq. (l24l >. so that expanded 
in terms of k gives the Hamiltonian 


H(X eq + + V(X eq , p 

2/le dp 






X=x„ 


+ o(/). 


(B18) 


Again, after including zero-point energy corrections, the error 
in the BOA ground state energy is 0(k 4 ) ~ m e /M, as shown 
in Fig. EJ 
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